Vocal interactions in Spix’s disc-winged bats (Thyroptera tricolor)
Author
Marcelo Araya-Salas
Published
September 20, 2023
Data analysis for the paper:
G Chaverri, M Sagot, JL Stynoski, M Araya-Salas, Y Araya-Ajoy, M Nagy, M Knörnschild, S Chaves-Ramírez, N Rose, M Sánchez-Chavarría, Y Jiménez-Torres, D Ulloa-Sanabria, H Solís-Hernández, GG Carter. In review. Contact-calling rates within groups of disc-winged bats vary by caller but not by the receiver’s identity, kinship, or association. Philosophical Transactions of the Royal Society B.
# set ggplot themetheme_set(theme_bw(base_size =14))###### Create functions to estimate confindence intervals###### function to get mean and 95% CI of values x via###### bootstrappingboot_ci <-function(x, perms =5000, bca = F) { get_mean <-function(x, d) {return(mean(x[d])) } x <-as.vector(na.omit(x)) mean <-mean(x)if (bca) { boot <-boot.ci(boot(data = x, statistic = get_mean, R = perms,parallel ="multicore", ncpus =4), type ="bca") low <- boot$bca[1, 4] high <- boot$bca[1, 5] } else { boot <-boot.ci(boot(data = x, statistic = get_mean, R = perms,parallel ="multicore", ncpus =4), type ="perc") low <- boot$perc[1, 4] high <- boot$perc[1, 5] }c(low = low, mean = mean, high = high, N =round(length(x)))}# function to get mean and 95% CI via bootstrapping of values y# within grouping variable xboot_ci2 <-function(d = d, y = d$y, x = d$x, perms =5000, bca = F) { df <-data.frame(effect =unique(x)) df$low <-NA df$mean <-NA df$high <-NA df$n.obs <-NAfor (i in1:nrow(df)) { ys <- y[which(x == df$effect[i])]if (length(ys) >1&var(ys) >0) { b <-boot_ci(y[which(x == df$effect[i])], perms = perms,bca = bca) df$low[i] <- b[1] df$mean[i] <- b[2] df$high[i] <- b[3] df$n.obs[i] <- b[4] } else { df$low[i] <-min(ys) df$mean[i] <-mean(ys) df$high[i] <-max(ys) df$n.obs[i] <-length(ys) } } df}# create function to convert matrix to data-framematrix_to_df <-function(m1) {data.frame(dyad =paste(rownames(m1)[col(m1)], colnames(m1)[row(m1)],sep ="_"), value =c(t(m1)), stringsAsFactors =FALSE)}# create function to plot variable by kinship with means and 95%# CIsplot_kinship_effect <-function(d = d, y = d$sri, label ="label") {# plot association by kinshipset.seed(123) means <- d %>%mutate(y = y) %>%filter(!is.na(kinship)) %>%group_by(dyad) %>%summarize(kinship =mean(kinship), y =mean(y), udyad.sex =first(udyad.sex)) %>%mutate(kin =ifelse(kinship >0, "kin", "nonkin")) %>% dplyr::select(kin, y) %>%boot_ci2(y = .$y, x = .$kin, bca = T) %>%rename(kin = effect) %>%mutate(kin =factor(kin, levels =c("nonkin", "kin"))) points <- d %>%mutate(y = y) %>%filter(!is.na(kinship)) %>%group_by(dyad) %>%summarize(kinship =mean(kinship), y =mean(y), udyad.sex =first(udyad.sex)) %>%mutate(kin =ifelse(kinship >0, "kin", "nonkin")) %>%mutate(kin =factor(kin, levels =c("nonkin", "kin")))# plot means and 95% CI means %>%ggplot(aes(x = kin, y = mean, color = kin)) +geom_jitter(data = points,aes(y = y), size =2, alpha =0.5, width =0.1) +geom_point(position =position_nudge(x =0.25),size =3) +geom_errorbar(aes(ymin = low, ymax = high, width =0.1),position =position_nudge(x =0.25), size =1) +xlab("") +ylab(label) +scale_color_viridis_d(alpha =0.8, begin =0.4,end =0.9, direction =-1, option ="G") +coord_flip() +theme_bw() +theme(legend.position ="none")}
1 Visualize social networks
1.1 Prepare data
Code
# Using# https://dshizuka.github.io/networkanalysis/example_make_sparrow_network.html# Input databatdat21 =read.csv("./data/raw/associations_2021.csv", header =TRUE,sep =";", colClasses =c("numeric", "character", "Date", "character","character", "character", "character", "character", "character","character", "character", "character", "character"))batdat22 =read.csv("./data/raw/associations_2022.csv", header =TRUE,sep =";", colClasses =c("numeric", "character", "Date", "character","character", "character", "character", "character", "character","character", "character"))# Add attribute dataattrib_21 <-read.csv("./data/raw/supp_21.csv", header =TRUE, sep =";",colClasses =c("character", "character"))attrib_22 <-read.csv("./data/raw/supp_22.csv", header =TRUE, sep =";",colClasses =c("character", "character"))# Prepare association data 2021batcols =grep("Bat", colnames(batdat21))# batcols batdat21[,batcols] gather(batdat21[,batcols])bat.ids =unique(gather(batdat21[, batcols])$value)# bat.idsbat.ids = bat.ids[!is.na(bat.ids)]# bat.idsbat.ids_df_21 <-data.frame(bat.ids)# csv_path <- 'bats_21.csv' write.csv(bat.ids_df_21, csv_path,# row.names = FALSE)m1_21 =apply(batdat21[, batcols], 1, function(x) match(bat.ids, x))# m1_21m1_21[is.na(m1_21)] =0m1_21[m1_21 >0] =1# m1_21rownames(m1_21) = bat.ids #rows are bat idscolnames(m1_21) =paste("group", 1:ncol(m1_21), sep ="_") #columns are group IDs (just 'group_#')# m1_21# rowSums(m1_21) set.seed(2) png(file='./output/Histogram# captures 2021.png', width=600, height=400) par(oma=c(1,1,1,1))# # all sides have 3 lines of space par(mar=c(5,4,4,2) + 0.8)# hist(rowSums(m1_21), main='', xlab='Number of observations',# ylab='Frequency', cex.lab = 2.5, cex.axis = 1.5) dev.off()average_n_21 <-mean(rowSums(m1_21)) #average number of times a single bat was observed# average_n_21sd_n_21 <-sd(rowSums(m1_21)) #sd number of times a single bat was observed# sd_n_21average_gs_21 <-mean(colSums(m1_21)) #average group size# average_gs_21sd_gs_21 <-sd(colSums(m1_21)) #sd group size# sd_gs_21# Prepare association data 2022batcols =grep("Bat", colnames(batdat22))# batcols batdat22[,batcols] gather(batdat22[,batcols])bat.ids =unique(gather(batdat22[, batcols])$value)# bat.idsbat.ids = bat.ids[is.na(bat.ids) == F]# bat.idsbat.ids_df_22 <-data.frame(bat.ids)# csv_path = 'bats_22.csv' write.csv(bat.ids_df_22, csv_path,# row.names = FALSE)m1_22 =apply(batdat22[, batcols], 1, function(x) match(bat.ids, x))# m1_22m1_22[is.na(m1_22)] =0m1_22[m1_22 >0] =1# m1_22rownames(m1_22) = bat.ids #rows are bat idscolnames(m1_22) =paste("group", 1:ncol(m1_22), sep ="_") #columns are group IDs (just 'group_#')# m1_22# rowSums(m1_22) set.seed(2) png(file='./output/Histogram# captures 2022.png', width=600, height=400) par(oma=c(1,1,1,1))# # all sides have 3 lines of space par(mar=c(5,4,4,2) + 0.8)# hist (rowSums(m1_22), main='', xlab='Number of observations',# ylab='Frequency', cex.lab = 2.5, cex.axis = 1.5) dev.off()average_n_22 <-mean(rowSums(m1_22)) #average number of times a single bat was observed# average_n_22sd_n_22 <-sd(rowSums(m1_22)) #sd number of times a single bat was observed# sd_n_22average_gs_22 <-mean(colSums(m1_22)) #average group size# average_gs_22sd_gs_22 <-sd(colSums(m1_22)) #average group size# sd_gs_22
1.2 Plot networks
Code
#Create adjacency data for associations (2021)adj_21=get_network(t(m1_21), data_format="GBI", association_index ="SRI") # the adjacency matrix
Generating 76 x 76 matrix
Code
g_21=graph_from_adjacency_matrix(adj_21, "undirected", weighted=T) #the igraph objectV(g_21)$sex <-factor(attrib_21[match(V(g_21)$name, attrib_21$bat_id), "sex"])mismatched_names <-setdiff(V(g_21)$name, attrib_21$bat_id) #check for issues#Create adjacency data for associations (2022)adj_22=get_network(t(m1_22), data_format="GBI", association_index ="SRI") # the adjacency matrix
Generating 76 x 76 matrix
Code
g_22=graph_from_adjacency_matrix(adj_22, "undirected", weighted=T) #the igraph objectV(g_22)$sex <-factor(attrib_22[match(V(g_22)$name, attrib_22$bat_id), "sex"])mismatched_names <-setdiff(V(g_22)$name, attrib_22$bat_id) #check for issues#Select individuals for which more than 10 observations were made (not incorporated into adjacency data for associations, above)#m2_21=m1_21[which(rowSums(m1_21)>10),]#adj_21=get_network(t(m2_21), data_format="GBI","SRI")#write.csv(adj_21, "assoc_21.csv", row.names=TRUE)#rowSums(m2_21)#m2_22=m1_22[which(rowSums(m1_22)>10),]#adj_22=get_network(t(m2_22), data_format="GBI","SRI")#rowSums(m2_22)#write.csv(adj_22, "assoc_22.csv", row.names=TRUE)#rowSums(m2_22)#Create communities (2021)com_21=fastgreedy.community(g_21) #community detection method# com_21[[11]]# Assign community namesmapping <-c(11, # Old membership 1 should be corrected to 118, # Old membership 2 should be corrected to 81, # Old membership 3 should be corrected to 12, # Old membership 4 should be corrected to 26, # Old membership 5 should be corrected to 69, # Old membership 6 should be corrected to 94, # Old membership 7 should be corrected to 410, # Old membership 8 should be corrected to 107, # Old membership 9 should be corrected to 75, # Old membership 10 should be corrected to 53# Old membership 11 should be corrected to 3)corrected_memberships <- mapping[com_21$membership]com_21$membership <- corrected_membershipscommunity_assignments <- com_21$membershipcolor_mapping <-c(mako(n =2, alpha =0.7, begin =0.4, end =0.9, direction =-1), "gray")names(color_mapping) <-c("female", "male", "unknown")node_sex <-V(g_21)$sexnode_colors <-sapply(node_sex, function(sex) { color_mapping[sex]})V(g_21)$color <- node_colors#Prepare to save figure # png(file="./output/Association networks 2021.png",# width=600, height=600)par(mfrow =c(1, 2), xpd =TRUE)# Iterate through com_21 and assign community names to nodes in g_21set.seed(123)plot( g_21, vertex.color = node_colors, # Use the calculated node colorsvertex.size =10, # Adjust the size of the nodes as neededvertex.label = community_assignments, # Use the extracted community assignments as labelsedge.width=E(g_21)$weight*10)title(main ="2021", cex.main =2)legend("topleft", legend =names(color_mapping), fill = color_mapping, bty ="n", inset=c(0.66, 0.6), x.intersp =0.2, y.intersp =0.55, cex =1.1)#Create communities (2022) com_22=fastgreedy.community(g_22) #community detection method# Assign community namesmapping <-c(3, 4, 1, 9, 10, 11, 7, 5, 6, 2, 12,8)corrected_memberships <- mapping[com_22$membership]com_22$membership <- corrected_membershipscommunity_assignments <- com_22$membershipnode_sex <-V(g_22)$sexnode_colors <-sapply(node_sex, function(sex) { color_mapping[sex]})V(g_22)$color <- node_colors#Prepare to save figure # png(file="./output/Association networks 2022.png",width=600, height=600)# Iterate through com_21 and assign community names to nodes in g_21set.seed(123)plot( g_22, vertex.color = node_colors, # Use the calculated node colorsvertex.size =10, # Adjust the size of the nodes as neededvertex.label = community_assignments, # Use the extracted community assignments as labelsedge.width=E(g_22)$weight*10)title(main ="2022", cex.main =2)
1.3 Estimate modulatity
Modularity:
2021: 0.87
2022: 0.9
2 Calling rate repeatability
2.1 Negative binomial model
Code
# Data manipulation Joining the two data setsd1 <-read.csv("./data/raw/vocal_interactions_2021.csv", sep =";")d2 <-read.csv("./data/raw/vocal_interactions_2022.csv", sep =";")colnames(d2)[6] <-"Dyad"d <-rbind(d1, d2)## Adding observation level random effectd$ob <-1:nrow(d)# Stats models Model for inquirymod_inquiry <-glmer.nb(n_inquiry ~1+ (1| Bat_roosting) + (1| Bat_flying) + (1| Group), data = d)summary(mod_inquiry)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: Negative Binomial(2.7219) ( log )
Formula: n_inquiry ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) + (1 |
Group)
Data: d
AIC BIC logLik deviance df.resid
4596.7 4617.2 -2293.3 4586.7 446
Scaled residuals:
Min 1Q Median 3Q Max
-1.62529 -0.48778 0.07216 0.47606 2.82516
Random effects:
Groups Name Variance Std.Dev.
Bat_roosting (Intercept) 1.523e-11 3.903e-06
Bat_flying (Intercept) 2.348e-01 4.845e-01
Group (Intercept) 2.141e-02 1.463e-01
Number of obs: 451, groups: Bat_roosting, 74; Bat_flying, 73; Group, 14
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.12062 0.07845 52.53 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see ?isSingular
# Code by Gerald Carter, gcarter1640@gmail.com### clean and wrangle data# get kinship datak <-as.matrix(read.csv('./data/raw/KinshipMatrix.csv', sep=",", row.names =1))# check symmetrymean(k==t(k), na.rm=T)# get sex datasex1 <-read.csv('supp_21.csv', sep=";") sex2 <-read.csv('supp_22.csv', sep=";") # get association dataa21 <-read.csv('associations_2021.csv', sep=";")a22 <-read.csv('associations_2022.csv', sep=";")# get calling datai21 <-read.csv('vocal_interactions_2021.csv', sep=";")i22 <-read.csv('vocal_interactions_2022.csv', sep=";")# tidy sex datasex <-rbind(sex1,sex2) %>%mutate(bat=paste0('bat', bat_id)) %>%group_by(bat, sex) %>%summarize(n=n()) %>% dplyr::select(bat, sex) %>%ungroup()# tidy 2021 association dataassoc21 <- a21 %>%pivot_longer(Bat1:Bat10, names_to ='names', values_to ="bat") %>%filter(!is.na(bat)) %>%mutate(bat=paste0('bat', bat)) %>%mutate(group=paste0('group', Group)) %>%mutate(date=dmy(Date)) %>% dplyr::select(obs= data_id, group, date, bat)# tidy 2022 association dataassoc22 <- a22 %>%pivot_longer(Bat1:Bat8, names_to ='names', values_to ="bat") %>%filter(!is.na(bat)) %>%mutate(bat=paste0('bat', bat)) %>%mutate(group=paste0('group', Group)) %>%mutate(date=dmy(Date)) %>% dplyr::select(obs= data_id, group, date, bat)# combine association dataassoc <-rbind(assoc21, assoc22) %>%mutate(obs=paste(obs,group,date, sep="_")) %>% dplyr::select(bat, obs) %>%as.data.frame()# get number of observations per batobs_per_bat <- assoc %>%group_by(bat) %>%summarize(n.obs=n()) %>%arrange(desc(n.obs))# range is 1 to 51 timesrange(obs_per_bat$n.obs)# pick threshold of observations for including bats in analysis# if they have fewer observations than we make their association rates NAthreshold <-4# plot number of observations per batassoc %>%group_by(bat) %>%summarize(n.obs=n()) %>%ggplot(aes(x=n.obs))+geom_histogram(fill="light blue", color="black")+geom_vline(xintercept= threshold, color='red', size=1)+xlab("number of observations of bat")+ylab("count of bats")# get bats seen fewer than threshold numberlow.n.bats <- assoc %>%group_by(bat) %>%summarize(n.obs=n()) %>%filter(n.obs < threshold) %>%pull(bat)low.n.batslength(low.n.bats)# 14 bats were seen fewer than 4 times### get SRI from asnipe# get association rates as group-by-individualgbi <-get_group_by_individual(assoc, data_format="individuals")# get association network of SRIassoc.net<-get_network(association_data=gbi, data_format ="GBI", association_index ="SRI")# check symmetrymean(assoc.net==t(assoc.net))# get SRI values for every undirected pair as dataframeSRI <-matrix_to_df(assoc.net) # tidy kinshipkinship <- k %>%matrix_to_df() %>%separate(dyad, into=c("bat1", "bat2")) %>%mutate(bat2=str_remove(bat2, "X")) %>%mutate(bat1=paste0('bat', bat1), bat2=paste0('bat', bat2)) %>%filter(bat1!=bat2) %>%# label undirected pairmutate(dyad=ifelse(bat1<bat2, paste(bat1,bat2, sep="_"), paste(bat2,bat1, sep="_"))) %>%group_by(dyad) %>%summarize(kinship=first(value)) # count trialsnrow(i21)nrow(i22)# trials with missing datarbind(i21,i22) %>%as_tibble() %>%filter(is.na(n_response)) %>%nrow()#2# trials where neither bat calledrbind(i21,i22) %>%as_tibble() %>%filter(n_inquiry==0) %>%nrow()# 5# fix and tidy calling datacalling <-rbind(i21,i22) %>%mutate(Date=as.character(mdy(Date))) %>%mutate(group=paste(substr(Date, start=1, stop=4),Group, sep="_")) %>%# need to recreate these labels because messed up by excel in raw dataseparate(Dyad, into=c("bat_flying", "bat_roosting"), remove=F) %>%mutate(bat_flying=paste0('bat', bat_flying), bat_roosting=paste0('bat', bat_roosting)) %>%# label undirected dyadsmutate(dyad=ifelse(bat_flying<bat_roosting, paste(bat_flying, bat_roosting, sep="_"), paste(bat_roosting, bat_flying, sep="_"))) %>% dplyr::select(date= Date, group, bat_flying, bat_roosting, dyad, flight_time, n_inquiry, n_response)# combine all datad <-# start with calling data calling %>%# add SRI mutate(sri= SRI$value[match(.$dyad, SRI$dyad)]) %>%# add kinshipmutate(kinship= kinship$kinship[match(.$dyad, kinship$dyad)]) %>%# add sexes of flying and roosting batsmutate(flying.sex= sex$sex[match(.$bat_flying, sex$bat)]) %>%mutate(roosting.sex= sex$sex[match(.$bat_roosting, sex$bat)]) %>%# label sex combination for flying-->roosting dyadmutate(dyad.sex=paste(flying.sex, roosting.sex, sep="-->")) %>%# label sex combination for dyad (male, female, both)mutate(udyad.sex =case_when( flying.sex =='female'& roosting.sex =='female'~"female-female", flying.sex =='male'& roosting.sex =='male'~"male-male",TRUE~"mixed-sex")) %>%# replace zero sri (never previously seen together) with NA because association is not 0mutate(sri =if_else(sri==0, NA, sri)) %>%# if flying bat seen fewer than 4 times, then SRI is NAmutate(sri =if_else(bat_flying %in% low.n.bats, NA, sri)) %>%# if roosting bat seen fewer than 4 times, then SRI is NAmutate(sri =if_else(bat_roosting %in% low.n.bats, NA, sri))# inspect and fix cases where flying bats in more than 2 groupst <- d %>%group_by(bat_flying, date, group) %>%summarize(n=n()) %>%arrange(date) %>%group_by(bat_flying, group) %>%summarize(n=n()) %>%group_by(bat_flying) %>%summarize(n=n()) %>%filter(n>2) %>%pull(bat_flying)d %>% dplyr::select(date, group, bat_flying) %>%filter(bat_flying %in% t) %>%arrange(bat_flying) %>%as.data.frame()# some of these seem impossible and must be errors# fix impossible group labelsd$group[which(d$date=="2021-07-03"& d$bat_flying=="bat982126051278521")] <-"2021_9"d$group[which(d$date=="2021-07-03"& d$bat_flying=="bat982126058484300")] <-"2021_9"# fix cases where roosting bats in more than 2 groupst <- d %>%group_by(bat_roosting, date, group) %>%summarize(n=n()) %>%arrange(date) %>%group_by(bat_roosting, group) %>%summarize(n=n()) %>%group_by(bat_roosting) %>%summarize(n=n()) %>%filter(n>2) %>%pull(bat_roosting)d %>% dplyr::select(date, group, bat_roosting) %>%filter(bat_roosting %in% t) %>%arrange(bat_roosting) %>%as.data.frame()# fix group labelsd$group[which(d$date=="2021-07-03"& d$bat_roosting=="bat982126058484300")] <-"2021_9"# group 9 split into groups 9,1 and 9,2 during 2022# There are interesting movements between group 9 and 10 and between groups 9,1 and 9,2# For this analysis, I'm going to consider group 9,1 and 9,2 as the same groupd$group[which(d$group=="2022_9,1"| d$group=="2022_9,2")] <-"2021_9"# remove trials without inquiry or response callsd <- d %>%# 453 rowsas_tibble() %>%filter(! (is.na(n_inquiry) &is.na(n_response))) %>%# 451 rowsfilter(n_inquiry>0) %>%# 446 rowsmutate(year=substr(date, 1,4)) # add year# save data as csvwrite.csv(d, file =paste0('./processed/calling-association-kinship-data.csv'), row.names= F)
3.2 Explore data
Code
# get dataraw <-read.csv("./data/processed/calling-association-kinship-data.csv")# look at data head(raw) group is the year and social group dyad# is the undirected pair (flying and roosting bat) in# alphanumeric order sri is simple ratio index of association# flying.sex is sex of flying bat roosting.sex is sex of# roosting bat dyad.sex is the sex of the flying and roosting# bat udyad.sex is females, males, or mixed# add a few variables to the datad <- raw %>%# get log count of inquiry callsmutate(log_inquiry =log(n_inquiry)) %>%# rescale kinship so 1 unit is 0.5mutate(kinship2 = kinship *2) %>%# if the roosting bat leaves the roost (escapes) then flight# time is < 300 smutate(escape =as.numeric(flight_time <300)) %>%# convert flight time from seconds to minutes (for easier# interpretation)mutate(flight_time2 = flight_time/60)###### Describe the data-------------# how many trials? d %>% nrow() 446# # how many undirected pairs have vocal data? d %>% pull(dyad)# %>% unique() %>% length 139 undirected pairs# d %>% group_by(bat_flying, bat_roosting) %>% summarize(n=n())# 254 directed pairs# how many groups? n_distinct(d$group) 23# how many of these have association data d %>% group_by(dyad)# %>% summarize(sri= mean(sri)) %>% filter(!is.na(sri)) 133# pairs have association# how many of these have kinship data d %>% group_by(dyad) %>%# summarize(kinship= mean(kinship)) %>% filter(!is.na(kinship))# 82 pairs have kinship data# how many related undirected pairs d %>% group_by(dyad) %>%# summarize(kinship= mean(kinship)) %>% filter(kinship>0) 32 kin# pairs# how many unrelated undirected pairs d %>% group_by(dyad) %>%# summarize(kinship= mean(kinship)) %>% filter(kinship==0) 50# nonkin pairs# what is mean kinship in group?grp.kin <- d %>%group_by(dyad, group) %>%summarize(kinship =mean(kinship, na.rm = T)) %>%group_by(group) %>%summarize(kinship =mean(kinship, na.rm = T), n =n()) %>%as.data.frame()set.seed(123)# grp.kin %>% pull(kinship) %>% boot_ci(bca=T) 0.23, 95% CI =# [0.16, 0.31]# how many groups have kin? grp.kin %>% filter(kinship>0) 17# how many groups have no kin? grp.kin %>% filter(kinship==0) 4# plot distribution of kinshipd %>%ggplot(aes(x = kinship)) +facet_wrap(~udyad.sex, ncol =1, scales ="free_y") +geom_histogram(fill =viridis(10, alpha =0.6)[3], color ="black") +ggtitle("pairwise kinship")
Code
# # check categories d %>% filter(!is.na(kinship)) %>%# group_by(dyad) %>% summarize(kinship= mean(kinship)) %>%# group_by(kinship) %>% summarize(n=n())##### Effect of kinship on association-------------# plot distribution of assocd %>%ggplot(aes(x = sri)) +facet_wrap(~udyad.sex, ncol =1, scale ="free_y") +geom_histogram(fill =viridis(10, alpha =0.6)[3], color ="black") +xlab("association rate") +ggtitle("pairwise SRI values")
Code
# looks ok# what is mean and 95% CI of assoc among flying and roosting# bats? set.seed(123) d %>% group_by(dyad) %>% summarize(sri =# mean(sri)) %>% pull(sri) %>% boot_ci(bca=T) 0.51, 95% CI =# [0.47, 0.55] a bit low, expected from past work is 0.76# now only nonkinset.seed(123)k1 <- d %>%filter(kinship ==0) %>%group_by(dyad) %>%summarize(sri =mean(sri)) %>%pull(sri) %>%boot_ci(bca = T) %>%c(kinship =0)# k1 0.572, 95% CI = [0.500, 0.636]# now only kin set.seed(123) d %>% filter(kinship>0) %>%# group_by(dyad) %>% summarize(sri = mean(sri)) %>% pull(sri)# %>% boot_ci(bca=T) 0.686, 95% CI = [0.631, 0.741]# now only close kinset.seed(123)k2 <- d %>%filter(kinship ==0.5) %>%group_by(dyad) %>%summarize(sri =mean(sri)) %>%pull(sri) %>%boot_ci(bca = T) %>%c(kinship =0.5)# k2 0.687, 95% CI = [0.628, 0.750]# now only 0.25 kinset.seed(123)k3 <- d %>%filter(kinship ==0.25) %>%group_by(dyad) %>%summarize(sri =mean(sri)) %>%pull(sri) %>%boot_ci(bca = T) %>%c(kinship =0.25)# k3 0.684, 95% CI = [0.542, 0.783]# compile meansk <-rbind(k1, k2, k3) %>%data.frame() %>%mutate(kin =as.character(kinship)) %>%mutate(kin2 = kinship >0) %>%mutate(assoc = mean)# plot association by kinship(kin.dyads.plot <- d %>%filter(!is.na(kinship)) %>%group_by(dyad) %>%summarize(kinship =mean(kinship), assoc =mean(sri), udyad.sex =first(udyad.sex)) %>%mutate(kin2 = kinship >0) %>%mutate(kin =as.character(kinship)) %>%ggplot(aes(x = kin, y = assoc, group = kin, color = kin2)) +geom_jitter(size =2,width =0.1, height =0, aes(shape = udyad.sex)) +geom_point(data = k,position =position_nudge(x =0.25), size =3) +geom_errorbar(data = k,aes(ymin = low, ymax = high, width =0.1), position =position_nudge(x =0.25),size =1) +xlab("kinship") +ylab("association rate (simple ratio index)") +scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1,option ="G") +theme(legend.position ="none", legend.title =element_blank()))
Code
# save as PDF ggsave( 'kin_association.pdf', plot =# kin.dyads.plot, scale = 1, width = 3, height = 5, units =# 'in', dpi = 300)##### Effect of flight time on calling ----------### How does flight time and count of inquiry calls predict### count of response calls?# plot number of response calls by flight timet1 <- d %>%mutate(kinship = kinship >0.1) %>%ggplot(aes(x = flight_time, y = n_response)) +geom_point(size =2,alpha =0.7, aes(color = kinship, shape = kinship)) +geom_smooth(method ="glm.nb",color =mako(10)[9]) +scale_color_viridis_d(alpha =1, begin =0.4,end =0.9, direction =-1, option ="G") +xlab("flight time (seconds)") +ylab("response call count") +theme(legend.position ="none")# plot number of inquiry calls by flight timet2 <- d %>%mutate(kinship = kinship >0.1) %>%ggplot(aes(x = flight_time, y = n_inquiry)) +geom_point(size =2,alpha =0.7, aes(color = kinship, shape = kinship)) +geom_smooth(method ="glm.nb",color =mako(10)[9]) +scale_color_viridis_d(alpha =1, begin =0.4,end =0.9, direction =-1, option ="G") +xlab("flight time (seconds)") +ylab("inquiry call count") +theme(legend.position ="none")# plot togethert <- t1 + t2t
# try poisson model for effect of flight time on inquiry callingt <-glmer(n_inquiry ~ flight_time2 + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d, family = poisson)check_overdispersion(t)
# Overdispersion test
dispersion ratio = 6.288
Pearson's Chi-Squared = 2773.031
p-value = < 0.001
Code
# negative binomial model (NBM)t <-glmmTMB(n_inquiry ~ flight_time2 + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d, ziformula =~0, family = nbinom2)check_overdispersion(t)
# Overdispersion test
dispersion ratio = 0.786
Pearson's Chi-Squared = 345.749
p-value = 1
# for every minute of flight time, the inquiry call count# increases by a factor of 1.43 [1.36, 1.49]# try poisson model for effect of flight time on response# callingt <-glmer(n_response ~ flight_time2 + log_inquiry + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d, family = poisson)check_overdispersion(t)
# Overdispersion test
dispersion ratio = 63.477
Pearson's Chi-Squared = 27930.070
p-value = < 0.001
# for every minute of flight time, the response call count# increases by a factor of 1.58 [1.26, 1.95]##### Effect of inquiry calling on response calling------------# plot with log countsd %>%ggplot(aes(x = log_inquiry, y =log(n_response +1))) +geom_point(size =2,color =adjustcolor("black", alpha.f =0.4)) +geom_smooth(method ="lm",color =mako(10)[9]) +xlab("log inquiry call count") +ylab("log (x+1) response call count")
Code
# plot negative binomial curve plot all batsd2 <- d %>%mutate(label ="all pairs") %>%mutate(escape =as.logical(escape))(t1 <-ggplot(d2, aes(x = log_inquiry, y = n_response)) +facet_wrap(~label) +geom_point(size =2, alpha =0.5) +geom_smooth(method ="glm.nb",color =mako(10)[9]) +xlab("log of inquiry call count") +ylab("response call count")) +theme(legend.position ="none")
Code
# plot with kinshipd3 <- d %>%filter(!is.na(kinship)) %>%mutate(escape =as.logical(escape)) %>%mutate(kin =if_else(kinship >0.1, "kin", "nonkin")) %>%mutate(kin =factor(kin, levels =c("nonkin", "kin")))d2$kin <-"all pairs"shared_cols <-intersect(names(d3), names(d2))dall <-rbind(d3[, shared_cols], d2[, shared_cols])dall$kin <-factor(dall$kin, levels =c("all pairs", "nonkin", "kin"))t2 <-ggplot(dall, aes(x = log_inquiry, y = (n_response), color = kin,group = kin)) +facet_wrap(~kin) +geom_point(size =2, alpha =0.5,aes(shape = escape)) +geom_smooth(method ="glm.nb", color =mako(10)[9]) +xlab("log of inquiry call count") +ylab("response call count") +scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1,option ="G") +theme(legend.position ="none")t2
Code
# save as PDF ggsave( 'inquiry_response.pdf', plot =# inquiry.response, width = 14, height = 7, units = 'in', dpi =# 300)# for plots, get residuals from negative binomial mixed effect# model predicting number of responses by number of inquiry# callsfit <-glmmTMB(n_response ~ log_inquiry +offset(log(flight_time)) + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d,ziformula =~1, family = nbinom2)# get difference between observed response count and predicted# response count add this response score to data for plotting# this value represents response calling more than expected from# inquiry callsd$resid_response <- d$n_response -predict(fit, type ="response",newdata = d, allow.new.levels = T)# plot number of inquiry calls by kinship with partnerd %>%filter(!is.na(kinship)) %>%mutate(kin =as.character(kinship)) %>%ggplot(aes(x = kin, y = n_inquiry, group = kin)) +geom_violin(fill =NA,width =1) +scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9,direction =-1, option ="G") +geom_jitter(size =2, width =0.1,height =0, aes(color = udyad.sex)) +xlab("kinship") +ylab("number of inquiry calls") +theme(legend.position ="top", legend.title =element_blank())
Code
# plot number of response calls by kinship with partnerd %>%filter(!is.na(kinship)) %>%mutate(kin =as.character(kinship)) %>%ggplot(aes(x = kin, y = n_response, group = kin)) +geom_violin(fill =NA,width =1) +geom_jitter(size =2, width =0.1, height =0, aes(color = udyad.sex)) +scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1,option ="G") +xlab("kinship") +ylab("number of response calls") +theme(legend.position ="top", legend.title =element_blank())
# plot effect of kinship on SRI(p1 <-plot_kinship_effect(d, d$sri, label ="association (SRI)"))
Code
# plot effect of kinship on inquiry calling(p2 <-plot_kinship_effect(d, d$n_inquiry/(d$flight_time/60), label ="inquiry calls per min of flight time"))
Code
# plot effect of kinship on response calling(p3 <-plot_kinship_effect(d, d$n_response/(d$flight_time/60), label ="response calls per min of flight time"))
Code
# plot effect of kinship on response calling(p4 <-plot_kinship_effect(d, d$resid_response, label ="residual response call variation"))
# fit zero-inflated NBMfit <-glmmTMB(n_response ~ kinship * sri + log_inquiry +offset(log(flight_time)) + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d,ziformula =~1, family = nbinom2)print("AIC")
[1] "AIC"
Code
AIC(fit) #2649
[1] 2648.867
Code
print("BIC")
[1] "BIC"
Code
BIC(fit) # 2687
[1] 2686.55
Code
summary(fit)
Family: nbinom2 ( log )
Formula:
n_response ~ kinship * sri + log_inquiry + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation: ~1
Data: d
AIC BIC logLik deviance df.resid
2648.9 2686.6 -1314.4 2628.9 310
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 3.792e-08 0.0001947
bat_roosting (Intercept) 2.296e-01 0.4791231
group (Intercept) 1.678e-01 0.4096551
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 0.707
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.2182 0.6294 -1.935 0.0529 .
kinship -3.3360 1.9159 -1.741 0.0817 .
sri -0.7478 0.6424 -1.164 0.2443
log_inquiry 0.1897 0.1154 1.644 0.1002
kinship:sri 5.2163 2.7487 1.898 0.0577 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4362 0.1242 -3.512 0.000444 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# fit zero-inflated NBM without interactionfit <-glmmTMB(n_response ~ kinship + sri + log_inquiry +offset(log(flight_time)) + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d,ziformula =~1, family = nbinom2)print("AIC")
[1] "AIC"
Code
AIC(fit) #2650
[1] 2650.47
Code
print("BIC")
[1] "BIC"
Code
BIC(fit) # 2684
[1] 2684.385
Code
summary(fit)
Family: nbinom2 ( log )
Formula:
n_response ~ kinship + sri + log_inquiry + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation: ~1
Data: d
AIC BIC logLik deviance df.resid
2650.5 2684.4 -1316.2 2632.5 311
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 4.139e-08 0.0002035
bat_roosting (Intercept) 1.976e-01 0.4445021
group (Intercept) 1.687e-01 0.4107856
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 0.678
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.6725 0.5816 -2.876 0.00403 **
kinship 0.2071 0.5173 0.400 0.68894
sri -0.2454 0.5926 -0.414 0.67882
log_inquiry 0.2338 0.1138 2.054 0.03995 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4456 0.1256 -3.548 0.000388 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# get predicted valuesd$n_response_predicted <-predict(fit, type ="response", newdata = d,allow.new.levels = T)# plot model performanced %>%filter(!is.na(kinship)) %>%mutate(kinship = kinship >0) %>%ggplot(aes(x = n_response, y = n_response_predicted)) +geom_point(size =2,aes(color = kinship)) +geom_smooth(method ="lm", color =mako(10)[9]) +xlab("observed count of response calls") +ylab("predicted count of response calls") +scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1,option ="G") +theme(legend.position ="top")
Code
# get relationship between predicted and observedsummary(lm(n_response_predicted ~ n_response, data = d)) # r-squared= 0.46
Call:
lm(formula = n_response_predicted ~ n_response, data = d)
Residuals:
Min 1Q Median 3Q Max
-105.210 -18.559 -1.269 18.411 109.081
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 58.98540 2.02304 29.16 <2e-16 ***
n_response 0.18176 0.01096 16.58 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 31.09 on 318 degrees of freedom
(126 observations deleted due to missingness)
Multiple R-squared: 0.4638, Adjusted R-squared: 0.4621
F-statistic: 275 on 1 and 318 DF, p-value: < 2.2e-16
Code
# get model coefficients with 95% CIscoefs <-tidy(fit, conf.int =TRUE, exponentiate = F, effects ="fixed",conf.method ="profile") %>%filter(component =="cond") %>% dplyr::select(-effect, -component) %>%mutate(type ="full model")# coefs# save model results write.csv(coefs,# file='response_model_results.csv')# effect of SRI with kin onlyfitk <-glmmTMB(n_response ~ sri + log_inquiry +offset(log(flight_time)) + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d %>%filter(kinship >0), ziformula =~1, family = nbinom2)coefs.k <-tidy(fitk, conf.int =TRUE, exponentiate = F, effects ="fixed",conf.method ="profile") %>%filter(component =="cond") %>% dplyr::select(-effect, -component) %>%mutate(type ="kin only")# coefs.k# effect of SRI with nonkin onlyfitn <-glmmTMB(n_response ~ sri + log_inquiry +offset(log(flight_time)) + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d %>%filter(kinship ==0), ziformula =~1, family = nbinom2)# get 95% CIcoefs.n <-tidy(fitn, conf.int =TRUE, exponentiate = F, effects ="fixed",conf.method ="profile") %>%filter(component =="cond") %>% dplyr::select(-effect, -component) %>%mutate(type ="nonkin only")# coefs.n# plot model resultstheme_set(theme_bw(base_size =12))plot1 <- coefs %>%filter(term !="(Intercept)") %>%mutate(term =case_when(term =="sri"~"association", term =="kinship"~"kinship", term =="log_inquiry"~"log count of inquiry calls", term =="kinship:sri"~"association x kinship interaction")) %>%mutate(term =factor(term, levels =c("log count of inquiry calls","association x kinship interaction", "association", "kinship"))) %>%mutate(term =fct_rev(term)) %>%ggplot(aes(x = estimate, y = term)) +geom_point(size =2) +geom_errorbarh(aes(xmin = conf.low,xmax = conf.high, height =0.2), size =1) +geom_vline(xintercept =0,linetype ="dashed") +ylab("") +xlab("Effect size") +coord_cartesian(xlim =c(-1.5,1.5)) +theme(axis.text =element_text(size =12), strip.text =element_text(size =12))plot1
Code
# save as PDF ggsave( 'response_model_results.pdf', plot =# plot1, scale = 1, width = 5, height = 2, units = 'in', dpi =# 300)# check that sri and kinship do not predict response calls as# single predictors fit zero-inflated negative binomial model# with only srit <-glmmTMB(n_response ~ sri + log_inquiry +offset(log(flight_time)) + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d,ziformula =~1, family = nbinom2)summary(t)
Family: nbinom2 ( log )
Formula:
n_response ~ sri + log_inquiry + offset(log(flight_time)) + (1 |
bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation: ~1
Data: d
AIC BIC logLik deviance df.resid
3579.0 3611.6 -1781.5 3563.0 428
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 2.229e-08 0.0001493
bat_roosting (Intercept) 6.142e-01 0.7836871
group (Intercept) 5.587e-02 0.2363637
Number of obs: 436, groups: bat_flying, 72; bat_roosting, 73; group, 23
Dispersion parameter for nbinom2 family (): 0.646
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.1155 0.5679 -3.725 0.000195 ***
sri -0.1597 0.4814 -0.332 0.740114
log_inquiry 0.3257 0.1102 2.955 0.003122 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4676 0.1155 -4.048 5.17e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# fit zero-inflated negative binomial model with only kinshipt <-glmmTMB(n_response ~ kinship + log_inquiry +offset(log(flight_time)) + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d,ziformula =~1, family = nbinom2)summary(t)
######## Effect of kinship and association on inquiry######## calling---------# poisson model is overdispersedt <-glmer(n_inquiry ~ kinship * sri +log(n_response +1) +offset(log(flight_time)) + (1| bat_flying) + (1| bat_roosting) + (1| group), family = poisson,data = d)check_overdispersion(t)
# Overdispersion test
dispersion ratio = 6.731
Pearson's Chi-Squared = 2100.164
p-value = < 0.001
Code
# fit negative binomial model (NBM) with log responsest <-glmmTMB(n_inquiry ~ kinship * sri +log(n_response +1) +offset(log(flight_time)) + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d,ziformula =~0, family = nbinom2)check_overdispersion(t)
# Overdispersion test
dispersion ratio = 0.820
Pearson's Chi-Squared = 255.132
p-value = 0.991
Code
AIC(t) # 3130
[1] 3129.521
Code
BIC(t) # 3163
[1] 3163.436
Code
# fit NBM with responsest <-glmmTMB(n_inquiry ~ kinship * sri + n_response +offset(log(flight_time)) + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d,ziformula =~0, family = nbinom2)check_overdispersion(t)
# Overdispersion test
dispersion ratio = 0.813
Pearson's Chi-Squared = 252.950
p-value = 0.993
Code
AIC(t) #3130
[1] 3130.339
Code
BIC(t) # 3164
[1] 3164.254
Code
summary(t)
Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship * sri + n_response + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
3130.3 3164.3 -1556.2 3112.3 311
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.833e-01 4.282e-01
bat_roosting (Intercept) 2.748e-09 5.242e-05
group (Intercept) 3.400e-02 1.844e-01
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 4.68
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3629986 0.1517410 -8.982 <2e-16 ***
kinship -0.5338331 0.6126168 -0.871 0.384
sri -0.1397377 0.1970890 -0.709 0.478
n_response 0.0001023 0.0002075 0.493 0.622
kinship:sri 0.4310938 0.9048669 0.476 0.634
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# fit NBM with responses without interactionfit2 <-glmmTMB(n_inquiry ~ kinship + sri + n_response +offset(log(flight_time)) + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d,ziformula =~0, family = nbinom2)AIC(fit2) #3129
[1] 3128.567
Code
BIC(fit2) # 3159
[1] 3158.713
Code
summary(fit2)
Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship + sri + n_response + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
3128.6 3158.7 -1556.3 3112.6 312
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.861e-01 0.4314410
bat_roosting (Intercept) 2.571e-09 0.0000507
group (Intercept) 3.441e-02 0.1855110
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 4.69
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3851696 0.1446519 -9.576 <2e-16 ***
kinship -0.2524912 0.1637266 -1.542 0.123
sri -0.1002280 0.1789141 -0.560 0.575
n_response 0.0001123 0.0002065 0.544 0.587
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# get predicted valuesd$n_inquiry_predicted <-predict(fit2, type ="response", newdata = d,allow.new.levels = T)# plot model performanced %>%filter(!is.na(kinship)) %>%mutate(kinship = kinship >0) %>%ggplot(aes(x = n_inquiry, y = n_inquiry_predicted)) +geom_point(size =2,aes(color = kinship)) +geom_smooth(method ="lm", color =mako(10)[9]) +xlab("observed count of inquiry calls") +ylab("predicted count of inquiry calls") +scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1,option ="G") +theme(legend.position ="top")
Code
# get relationship between predicted and observedsummary(lm(n_inquiry_predicted ~ n_inquiry, data = d)) # r-squared= 0.71
Call:
lm(formula = n_inquiry_predicted ~ n_inquiry, data = d)
Residuals:
Min 1Q Median 3Q Max
-92.258 -12.967 -0.546 11.381 87.999
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.40700 2.06189 9.412 <2e-16 ***
n_inquiry 0.71686 0.02554 28.063 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19.59 on 318 degrees of freedom
(126 observations deleted due to missingness)
Multiple R-squared: 0.7124, Adjusted R-squared: 0.7114
F-statistic: 787.5 on 1 and 318 DF, p-value: < 2.2e-16
Code
# get model coefficientssummary(fit2)
Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship + sri + n_response + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
3128.6 3158.7 -1556.3 3112.6 312
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.861e-01 0.4314410
bat_roosting (Intercept) 2.571e-09 0.0000507
group (Intercept) 3.441e-02 0.1855110
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 4.69
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3851696 0.1446519 -9.576 <2e-16 ***
kinship -0.2524912 0.1637266 -1.542 0.123
sri -0.1002280 0.1789141 -0.560 0.575
n_response 0.0001123 0.0002065 0.544 0.587
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# save as PDF ggsave( 'inquiry_model_results.pdf', plot = plot1,# scale = 1, width = 7, height = 3.5, units = 'in', dpi = 300)# get relative amount of variance explained by random intercepttidy(fit2, conf.int = F, exponentiate = F) %>%filter(effect =="ran_pars") %>%mutate(variance = estimate^2) %>% dplyr::select(effect, group, variance) %>%group_by(effect) %>%mutate(sum =sum(variance)) %>%mutate(prop = variance/sum)
# A tibble: 3 × 5
# Groups: effect [1]
effect group variance sum prop
<chr> <chr> <dbl> <dbl> <dbl>
1 ran_pars bat_flying 0.186 0.221 0.844
2 ran_pars bat_roosting 0.00000000257 0.221 0.0000000117
3 ran_pars group 0.0344 0.221 0.156
Code
# confirm no effect with single predictors associationt <-glmmTMB(n_inquiry ~ sri + n_response +offset(log(flight_time)) + (1| bat_flying) + (1| bat_roosting) + (1| group), data = d,ziformula =~0, family = nbinom2)summary(t)
Family: nbinom2 ( log )
Formula:
n_inquiry ~ sri + n_response + offset(log(flight_time)) + (1 |
bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
4257.0 4285.5 -2121.5 4243.0 429
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 0.23098 0.48060
bat_roosting (Intercept) 0.00448 0.06694
group (Intercept) 0.01246 0.11164
Number of obs: 436, groups: bat_flying, 72; bat_roosting, 73; group, 23
Dispersion parameter for nbinom2 family (): 5.1
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3782473 0.1046922 -13.165 <2e-16 ***
sri -0.1730846 0.1317072 -1.314 0.1888
n_response 0.0003027 0.0001691 1.790 0.0734 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---title: Data analysissubtitle: Vocal interactions in Spix's disc-winged bats (*Thyroptera tricolor*)author: Marcelo Araya-Salasdate: "`r Sys.Date()`"toc: truetoc-depth: 3toc-location: leftnumber-sections: truehighlight-style: pygmentsformat: html: code-fold: true code-tools: true code-copy: true embed-resources: trueeditor_options: chunk_output_type: console---<!-- skyblue box --><divclass="alert alert-info">Data analysis for the paper:G Chaverri, M Sagot, JL Stynoski, M Araya-Salas, Y Araya-Ajoy, M Nagy, M Knörnschild, S Chaves-Ramírez, N Rose, M Sánchez-Chavarría, Y Jiménez-Torres, D Ulloa-Sanabria, H Solís-Hernández, GG Carter. In review. **Contact-calling rates within groups of disc-winged bats vary by caller but not by the receiver’s identity, kinship, or association**. Philosophical Transactions of the Royal Society B. </div> <!-- this code add line numbers to code blocks --><!-- only works when code folding is not used in yaml (code_folding: show) -->```{=html}<style>body { counter-reset: source-line 0; }pre.numberSource code { counter-reset: none; }</style>``````{r set root directory, echo = FALSE}# set working directory as project directory or one directory above,rootdir <-try(rprojroot::find_rstudio_root_file(), silent =TRUE)if (is(rootdir, "try-error")) rootdir <-".."knitr::opts_knit$set(root.dir = rootdir)``````{r add link to github repo, echo = FALSE, results='asis'}# print link to github repo if anyif (file.exists("./.git/config")){ config <-readLines("./.git/config") url <-grep("url", config, value =TRUE) url <-gsub("\\turl = |.git$", "", url)[2]cat("\nSource code and data found at [", url, "](", url, ")", sep ="") }``````{r setup style, echo = FALSE, message = FALSE, warning=FALSE}# options to customize chunk outputsknitr::opts_chunk$set(class.source ="numberLines lineAnchors", # for code line numberstidy.opts =list(width.cutoff =65), tidy =TRUE,message =FALSE,warning =FALSE )knitr::opts_knit$set(root.dir ="..")```# Load packages {.unnumbered .unlisted}```{r load packages}# load function from sketchysource("https://raw.githubusercontent.com/maRce10/sketchy/main/R/load_packages.R")# install/ load packagessketchy::load_packages(packages =c("igraph", "bipartite", "asnipe", "assortnet", "ggplot2", "ggmap", "ecodist", "igraphdata", "statnet", "RColorBrewer", "tidyverse", "geosphere", "mapproj", "sna", "stringr", "lme4", "glmmTMB", "broom.mixed", "boot", "patchwork", "performance", "MASS", "viridis"))```# Functions and global parameters {.unnumbered .unlisted}```{r}# set ggplot themetheme_set(theme_bw(base_size =14))###### Create functions to estimate confindence intervals# function to get mean and 95% CI of values x via bootstrappingboot_ci <-function(x, perms=5000, bca=F) { get_mean <-function(x, d) {return(mean(x[d])) } x <-as.vector(na.omit(x)) mean <-mean(x)if(bca){ boot <-boot.ci(boot(data=x, statistic=get_mean, R=perms, parallel ="multicore", ncpus =4), type="bca") low <- boot$bca[1,4] high <- boot$bca[1,5] }else{ boot <-boot.ci(boot(data=x, statistic=get_mean, R=perms, parallel ="multicore", ncpus =4), type="perc") low <- boot$perc[1,4] high <- boot$perc[1,5] }c(low=low,mean=mean,high=high, N=round(length(x)))}# function to get mean and 95% CI via bootstrapping of values y within grouping variable xboot_ci2 <-function(d=d, y=d$y, x=d$x, perms=5000, bca=F){ df <-data.frame(effect=unique(x)) df$low <-NA df$mean <-NA df$high <-NA df$n.obs <-NAfor (i in1:nrow(df)) { ys <- y[which(x==df$effect[i])]if (length(ys)>1&var(ys)>0 ){ b <-boot_ci(y[which(x==df$effect[i])], perms=perms, bca=bca) df$low[i] <- b[1] df$mean[i] <- b[2] df$high[i] <- b[3] df$n.obs[i] <- b[4] }else{ df$low[i] <-min(ys) df$mean[i] <-mean(ys) df$high[i] <-max(ys) df$n.obs[i] <-length(ys) } } df}# create function to convert matrix to data-framematrix_to_df <-function(m1){data.frame(dyad =paste(rownames(m1)[col(m1)], colnames(m1)[row(m1)], sep="_"),value =c(t(m1)), stringsAsFactors =FALSE)}# create function to plot variable by kinship with means and 95% CIsplot_kinship_effect <-function(d=d, y= d$sri, label='label') {# plot association by kinshipset.seed(123) means <- d %>%mutate(y=y) %>%filter(!is.na(kinship)) %>%group_by(dyad) %>%summarize(kinship=mean(kinship), y=mean(y), udyad.sex=first(udyad.sex)) %>%mutate(kin=ifelse(kinship>0, "kin", "nonkin")) %>% dplyr::select(kin, y) %>%boot_ci2(y= .$y, x=.$kin, bca=T) %>%rename(kin= effect) %>%mutate(kin=factor(kin, levels=c("nonkin", "kin"))) points <- d %>%mutate(y=y) %>%filter(!is.na(kinship)) %>%group_by(dyad) %>%summarize(kinship=mean(kinship), y=mean(y), udyad.sex=first(udyad.sex)) %>%mutate(kin=ifelse(kinship>0, "kin", "nonkin")) %>%mutate(kin=factor(kin, levels=c("nonkin", "kin")))# plot means and 95% CI means %>%ggplot(aes(x=kin, y=mean, color= kin))+geom_jitter(data= points, aes(y= y), size=2, alpha=0.5, width=0.1)+geom_point(position =position_nudge(x =0.25), size=3)+geom_errorbar(aes(ymin=low, ymax=high, width=.1), position =position_nudge(x =0.25), size=1)+xlab("")+ylab(label) +scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1, option ="G") +coord_flip()+theme_bw()+theme(legend.position="none")}```# Visualize social networks## Prepare data```{r}#Using https://dshizuka.github.io/networkanalysis/example_make_sparrow_network.html#Input databatdat21=read.csv('./data/raw/associations_2021.csv', header =TRUE, sep =";", colClasses =c("numeric", "character", "Date", "character", "character", "character","character","character","character","character","character","character","character"))batdat22=read.csv('./data/raw/associations_2022.csv', header =TRUE, sep =";", colClasses =c("numeric", "character", "Date", "character", "character", "character","character","character","character","character","character"))#Add attribute dataattrib_21 <-read.csv('./data/raw/supp_21.csv', header =TRUE, sep =";", colClasses =c("character", "character"))attrib_22 <-read.csv('./data/raw/supp_22.csv', header =TRUE, sep =";", colClasses =c("character", "character"))#Prepare association data 2021batcols=grep("Bat",colnames(batdat21))# batcols# batdat21[,batcols]# gather(batdat21[,batcols])bat.ids=unique(gather(batdat21[,batcols])$value)# bat.idsbat.ids=bat.ids[!is.na(bat.ids)]# bat.idsbat.ids_df_21 <-data.frame(bat.ids)# csv_path <- 'bats_21.csv'# write.csv(bat.ids_df_21, csv_path, row.names = FALSE)m1_21=apply(batdat21[,batcols], 1, function(x) match(bat.ids,x))# m1_21m1_21[is.na(m1_21)]=0m1_21[m1_21>0]=1# m1_21rownames(m1_21)=bat.ids #rows are bat idscolnames(m1_21)=paste('group', 1:ncol(m1_21), sep="_") #columns are group IDs (just "group_#")# m1_21# rowSums(m1_21)# set.seed(2)# png(file="./output/Histogram captures 2021.png", width=600, height=400)# par(oma=c(1,1,1,1)) # all sides have 3 lines of space# par(mar=c(5,4,4,2) + 0.8)# hist(rowSums(m1_21), main="", xlab="Number of observations", ylab="Frequency", cex.lab = 2.5, cex.axis = 1.5)# dev.off()average_n_21 <-mean(rowSums(m1_21)) #average number of times a single bat was observed# average_n_21sd_n_21 <-sd(rowSums(m1_21)) #sd number of times a single bat was observed# sd_n_21average_gs_21 <-mean(colSums(m1_21)) #average group size# average_gs_21sd_gs_21 <-sd(colSums(m1_21)) #sd group size# sd_gs_21#Prepare association data 2022batcols=grep("Bat",colnames(batdat22))# batcols# batdat22[,batcols]# gather(batdat22[,batcols])bat.ids=unique(gather(batdat22[,batcols])$value)# bat.idsbat.ids=bat.ids[is.na(bat.ids)==F]# bat.idsbat.ids_df_22 <-data.frame(bat.ids)# csv_path = 'bats_22.csv'# write.csv(bat.ids_df_22, csv_path, row.names = FALSE)m1_22=apply(batdat22[,batcols], 1, function(x) match(bat.ids,x))# m1_22m1_22[is.na(m1_22)]=0m1_22[m1_22>0]=1# m1_22rownames(m1_22)=bat.ids #rows are bat idscolnames(m1_22)=paste('group', 1:ncol(m1_22), sep="_") #columns are group IDs (just "group_#")# m1_22# rowSums(m1_22)# set.seed(2)# png(file="./output/Histogram captures 2022.png",# width=600, height=400)# par(oma=c(1,1,1,1)) # all sides have 3 lines of space# par(mar=c(5,4,4,2) + 0.8)# hist (rowSums(m1_22), main="", xlab="Number of observations", ylab="Frequency", cex.lab = 2.5, cex.axis = 1.5)# dev.off()average_n_22 <-mean(rowSums(m1_22)) #average number of times a single bat was observed# average_n_22sd_n_22 <-sd(rowSums(m1_22)) #sd number of times a single bat was observed# sd_n_22average_gs_22 <-mean(colSums(m1_22)) #average group size# average_gs_22sd_gs_22 <-sd(colSums(m1_22)) #average group size# sd_gs_22```## Plot networks```{r plot networks, message=FALSE, out.width= "100%", out.height="100%", fig.width=12, fig.height=12}#Create adjacency data for associations (2021)adj_21=get_network(t(m1_21), data_format="GBI", association_index ="SRI") # the adjacency matrixg_21=graph_from_adjacency_matrix(adj_21, "undirected", weighted=T) #the igraph objectV(g_21)$sex <-factor(attrib_21[match(V(g_21)$name, attrib_21$bat_id), "sex"])mismatched_names <-setdiff(V(g_21)$name, attrib_21$bat_id) #check for issues#Create adjacency data for associations (2022)adj_22=get_network(t(m1_22), data_format="GBI", association_index ="SRI") # the adjacency matrixg_22=graph_from_adjacency_matrix(adj_22, "undirected", weighted=T) #the igraph objectV(g_22)$sex <-factor(attrib_22[match(V(g_22)$name, attrib_22$bat_id), "sex"])mismatched_names <-setdiff(V(g_22)$name, attrib_22$bat_id) #check for issues#Select individuals for which more than 10 observations were made (not incorporated into adjacency data for associations, above)#m2_21=m1_21[which(rowSums(m1_21)>10),]#adj_21=get_network(t(m2_21), data_format="GBI","SRI")#write.csv(adj_21, "assoc_21.csv", row.names=TRUE)#rowSums(m2_21)#m2_22=m1_22[which(rowSums(m1_22)>10),]#adj_22=get_network(t(m2_22), data_format="GBI","SRI")#rowSums(m2_22)#write.csv(adj_22, "assoc_22.csv", row.names=TRUE)#rowSums(m2_22)#Create communities (2021)com_21=fastgreedy.community(g_21) #community detection method# com_21[[11]]# Assign community namesmapping <-c(11, # Old membership 1 should be corrected to 118, # Old membership 2 should be corrected to 81, # Old membership 3 should be corrected to 12, # Old membership 4 should be corrected to 26, # Old membership 5 should be corrected to 69, # Old membership 6 should be corrected to 94, # Old membership 7 should be corrected to 410, # Old membership 8 should be corrected to 107, # Old membership 9 should be corrected to 75, # Old membership 10 should be corrected to 53# Old membership 11 should be corrected to 3)corrected_memberships <- mapping[com_21$membership]com_21$membership <- corrected_membershipscommunity_assignments <- com_21$membershipcolor_mapping <-c(mako(n =2, alpha =0.7, begin =0.4, end =0.9, direction =-1), "gray")names(color_mapping) <-c("female", "male", "unknown")node_sex <-V(g_21)$sexnode_colors <-sapply(node_sex, function(sex) { color_mapping[sex]})V(g_21)$color <- node_colors#Prepare to save figure # png(file="./output/Association networks 2021.png",# width=600, height=600)par(mfrow =c(1, 2), xpd =TRUE)# Iterate through com_21 and assign community names to nodes in g_21set.seed(123)plot( g_21, vertex.color = node_colors, # Use the calculated node colorsvertex.size =10, # Adjust the size of the nodes as neededvertex.label = community_assignments, # Use the extracted community assignments as labelsedge.width=E(g_21)$weight*10)title(main ="2021", cex.main =2)legend("topleft", legend =names(color_mapping), fill = color_mapping, bty ="n", inset=c(0.66, 0.6), x.intersp =0.2, y.intersp =0.55, cex =1.1)#Create communities (2022) com_22=fastgreedy.community(g_22) #community detection method# Assign community namesmapping <-c(3, 4, 1, 9, 10, 11, 7, 5, 6, 2, 12,8)corrected_memberships <- mapping[com_22$membership]com_22$membership <- corrected_membershipscommunity_assignments <- com_22$membershipnode_sex <-V(g_22)$sexnode_colors <-sapply(node_sex, function(sex) { color_mapping[sex]})V(g_22)$color <- node_colors#Prepare to save figure # png(file="./output/Association networks 2022.png",width=600, height=600)# Iterate through com_21 and assign community names to nodes in g_21set.seed(123)plot( g_22, vertex.color = node_colors, # Use the calculated node colorsvertex.size =10, # Adjust the size of the nodes as neededvertex.label = community_assignments, # Use the extracted community assignments as labelsedge.width=E(g_22)$weight*10)title(main ="2022", cex.main =2)```## Estimate modulatityModularity:- 2021: `r round(modularity(com_21), 2)`- 2022: `r round(modularity(com_22), 2)`# Calling rate repeatability## Negative binomial model```{r}#Data manipulation## Joining the two data setsd1<-read.csv("./data/raw/vocal_interactions_2021.csv", sep=";")d2<-read.csv("./data/raw/vocal_interactions_2022.csv", sep=";")colnames(d2)[6]<-"Dyad"d<-rbind(d1,d2)## Adding observation level random effectd$ob<-1:nrow(d)#Stats models##Model for inquirymod_inquiry<-glmer.nb(n_inquiry~1+ (1|Bat_roosting) + (1|Bat_flying) + (1|Group), data=d)summary(mod_inquiry)###Extracting variance componentsmu_l<-fixef(mod_inquiry)[1]mu<-exp(mu_l)theta<-getME(mod_inquiry, "glmer.nb.theta")VBF_I_l<-VarCorr(mod_inquiry)$Bat_flying[1,1] #Variance Bats FlyingVBF_I<- (exp(VBF_I_l)-1)*exp(2*mu_l + VBF_I_l)VBR_I_l<-VarCorr(mod_inquiry)$Bat_roosting[1,1] #Variance Bats RoostingVBR_I<-(exp(VBR_I_l)-1)*exp(2*mu_l + VBR_I_l)VBG_I_l<-VarCorr(mod_inquiry)$Group[1,1] #Variance GRoupVBG_I<-(exp(VBG_I_l)-1)*exp(2*mu_l + VBG_I_l)NBV_I<-mu + mu^2/theta####Repeatabilty roosting individualprint("Repeatabilty roosting individual")VBR_I/(VBR_I+VBF_I+NBV_I + VBG_I)####Repeatabilty flying individualprint("Repeatabilty flying individual")VBF_I/(VBR_I+VBF_I+NBV_I)#Stats models##Model for inquiry2d2<-d[d$n_response>0,]mod_inquiry<-glmer.nb(n_inquiry~1+ (1|Bat_roosting)+ (1|Bat_flying) + (1|Group), data=d2)summary(mod_inquiry)###Extracting variance componentsmu_l<-fixef(mod_inquiry)[1]mu<-exp(mu_l)theta<-getME(mod_inquiry, "glmer.nb.theta")VBF_I_l<-VarCorr(mod_inquiry)$Bat_flying[1,1] #Variance Bats FlyingVBF_I<- (exp(VBF_I_l)-1)*exp(2*mu_l + VBF_I_l)VBR_I_l<-VarCorr(mod_inquiry)$Bat_roosting[1,1] #Variance Bats RoostingVBR_I<-(exp(VBR_I_l)-1)*exp(2*mu_l + VBR_I_l)VBG_I_l<-VarCorr(mod_inquiry)$Group[1,1] #Variance GRoupVBG_I<-(exp(VBG_I_l)-1)*exp(2*mu_l + VBG_I_l)NBV_I<-mu + mu^2/theta####Repeatabilty roosting individualprint("Repeatabilty roosting individual")VBR_I/(VBR_I+VBF_I+NBV_I)####Repeatabilty flying individualprint("Repeatabilty flying individual")VBF_I/(VBR_I+VBF_I+NBV_I)##Model for responsemod_response<-glmer.nb(n_response~1+ (1|Bat_roosting)+ (1|Bat_flying) + (1|Group), data=d)summary(mod_response)###Extracting variance componentsmu_l<-fixef(mod_response)[1] #+ fixef(mod_response)[2]*mean(d$n_inquiry, na.rm=TRUE)mu<-exp(mu_l)theta<-getME(mod_response, "glmer.nb.theta")VBF_R_l<-VarCorr(mod_response)$Bat_flying[1,1] #Variance Bats FlyingVBF_R<- (exp(VBF_R_l)-1)*exp(2*mu_l + VBF_R_l)VBR_R_l<-VarCorr(mod_response)$Bat_roosting[1,1] #Variance Bats RoostingVBR_R<-(exp(VBR_R_l)-1)*exp(2*mu_l + VBR_R_l)VBG_R_l<-VarCorr(mod_response)$Group[1,1] #Variance GRoupVBG_R<-(exp(VBG_R_l)-1)*exp(2*mu_l + VBG_R_l)NBV_R<-mu + mu^2/theta####Repeatabilty roosting individualprint("Repeatabilty roosting individual")VBR_R/(VBR_R+VBF_R+NBV_R +VBG_R)####Repeatabilty flying individualprint("Repeatabilty flying individual")VBF_R/(VBR_R+VBF_R+NBV_R+VBG_R)```## Poisson model```{r}# #Data manipulation# ## Joining the two data sets# d1<-read.csv("vocal_interactions_2021.csv", sep=";")# d2<-read.csv("vocal_interactions_2022.csv", sep=";")# # colnames(d2)[6]<-"Dyad"# d<-rbind(d1,d2)## Adding observation level random effectd$ob<-1:nrow(d)#Stats models##Model for inquiry with all observationsmod_inquiry<-glmer(n_inquiry~1+ (1|Bat_roosting)+ (1|Bat_flying) + (1|Group) + (1|ob), data=d, family="poisson")###Extracting variance componentsVBF_I<-VarCorr(mod_inquiry)$Bat_flying[1,1] #Variance Bats FlyingVBR_I<-VarCorr(mod_inquiry)$Bat_roosting[1,1] #Variance Bats RoostingVBG_I<-VarCorr(mod_inquiry)$Group[1,1] #Variance GroupVOD_I<-VarCorr(mod_inquiry)$ob[1,1] #Variance OverddispersionPV_I<-log(1/exp(fixef(mod_inquiry)[1]) +1) #Variance poisson processs####Repeatabilty roosting individualVBR_I/(VBR_I+VOD_I +VBF_I+PV_I + VBG_I)####Repeatabilty flying individualVBF_I/(VBR_I+VOD_I +VBF_I+PV_I+ VBG_I)##Model for inquiry only when someone respondedd2<-d[d$n_response>0,]mod_inquiry2<-glmer(n_inquiry~1+ (1|Bat_roosting)+ (1|Bat_flying) + (1|ob), data=d2, family="poisson")###Extracting variance componentsVBF_I<-VarCorr(mod_inquiry2)$Bat_flying[1,1] #Variance Bats FlyingVBR_I<-VarCorr(mod_inquiry2)$Bat_roosting[1,1] #Variance Bats RoostingVBG_I<-VarCorr(mod_inquiry)$Group[1,1]VOD_I<-VarCorr(mod_inquiry2)$ob[1,1] #Variance OverddispersionPV_I<-log(1/exp(fixef(mod_inquiry2)[1]) +1) #Variance poisson processs####Repeatabilty roosting individualprint("Repeatabilty roosting individual")VBR_I/(VBR_I+VOD_I +VBF_I+PV_I+ VBG_I)####Repeatabilty flying individualprint("Repeatabilty flying individual")VBF_I/(VBR_I+VOD_I +VBF_I+PV_I+ VBG_I)##Model for response mod_response<-glmer(n_response~1+ (1|Bat_roosting)+ (1|Bat_flying) + (1|Group) + (1|ob), data=d, family="poisson")summary(mod_response)###Extracting variance componentsVBF_R<-VarCorr(mod_response)$Bat_flying[1,1] #Variance Bats FlyingVBR_R<-VarCorr(mod_response)$Bat_roosting[1,1] #Variance Bats RoostingVBG_R<-VarCorr(mod_response)$Group[1,1]VOD_R<-VarCorr(mod_response)$ob[1,1] #Variance OverddispersionPV_R<-log(1/exp(fixef(mod_response)[1]) +1) #Variance poisson processs##Repeatabilty roosting individualprint("Repeatabilty roosting individual")VBR_R/(VBR_R + VOD_R + VBF_R + PV_R + VBG_R)##Repeatabilty flying individualprint("Repeatabilty flying individual")VBF_R/(VBR_R+ VOD_R + VBF_R + PV_R + VBG_R)```# Calling, association and kinship## Prepare data```{r, eval = FALSE}# Code by Gerald Carter, gcarter1640@gmail.com### clean and wrangle data# get kinship datak <-as.matrix(read.csv('./data/raw/KinshipMatrix.csv', sep=",", row.names =1))# check symmetrymean(k==t(k), na.rm=T)# get sex datasex1 <-read.csv('supp_21.csv', sep=";") sex2 <-read.csv('supp_22.csv', sep=";") # get association dataa21 <-read.csv('associations_2021.csv', sep=";")a22 <-read.csv('associations_2022.csv', sep=";")# get calling datai21 <-read.csv('vocal_interactions_2021.csv', sep=";")i22 <-read.csv('vocal_interactions_2022.csv', sep=";")# tidy sex datasex <-rbind(sex1,sex2) %>%mutate(bat=paste0('bat', bat_id)) %>%group_by(bat, sex) %>%summarize(n=n()) %>% dplyr::select(bat, sex) %>%ungroup()# tidy 2021 association dataassoc21 <- a21 %>%pivot_longer(Bat1:Bat10, names_to ='names', values_to ="bat") %>%filter(!is.na(bat)) %>%mutate(bat=paste0('bat', bat)) %>%mutate(group=paste0('group', Group)) %>%mutate(date=dmy(Date)) %>% dplyr::select(obs= data_id, group, date, bat)# tidy 2022 association dataassoc22 <- a22 %>%pivot_longer(Bat1:Bat8, names_to ='names', values_to ="bat") %>%filter(!is.na(bat)) %>%mutate(bat=paste0('bat', bat)) %>%mutate(group=paste0('group', Group)) %>%mutate(date=dmy(Date)) %>% dplyr::select(obs= data_id, group, date, bat)# combine association dataassoc <-rbind(assoc21, assoc22) %>%mutate(obs=paste(obs,group,date, sep="_")) %>% dplyr::select(bat, obs) %>%as.data.frame()# get number of observations per batobs_per_bat <- assoc %>%group_by(bat) %>%summarize(n.obs=n()) %>%arrange(desc(n.obs))# range is 1 to 51 timesrange(obs_per_bat$n.obs)# pick threshold of observations for including bats in analysis# if they have fewer observations than we make their association rates NAthreshold <-4# plot number of observations per batassoc %>%group_by(bat) %>%summarize(n.obs=n()) %>%ggplot(aes(x=n.obs))+geom_histogram(fill="light blue", color="black")+geom_vline(xintercept= threshold, color='red', size=1)+xlab("number of observations of bat")+ylab("count of bats")# get bats seen fewer than threshold numberlow.n.bats <- assoc %>%group_by(bat) %>%summarize(n.obs=n()) %>%filter(n.obs < threshold) %>%pull(bat)low.n.batslength(low.n.bats)# 14 bats were seen fewer than 4 times### get SRI from asnipe# get association rates as group-by-individualgbi <-get_group_by_individual(assoc, data_format="individuals")# get association network of SRIassoc.net<-get_network(association_data=gbi, data_format ="GBI", association_index ="SRI")# check symmetrymean(assoc.net==t(assoc.net))# get SRI values for every undirected pair as dataframeSRI <-matrix_to_df(assoc.net) # tidy kinshipkinship <- k %>%matrix_to_df() %>%separate(dyad, into=c("bat1", "bat2")) %>%mutate(bat2=str_remove(bat2, "X")) %>%mutate(bat1=paste0('bat', bat1), bat2=paste0('bat', bat2)) %>%filter(bat1!=bat2) %>%# label undirected pairmutate(dyad=ifelse(bat1<bat2, paste(bat1,bat2, sep="_"), paste(bat2,bat1, sep="_"))) %>%group_by(dyad) %>%summarize(kinship=first(value)) # count trialsnrow(i21)nrow(i22)# trials with missing datarbind(i21,i22) %>%as_tibble() %>%filter(is.na(n_response)) %>%nrow()#2# trials where neither bat calledrbind(i21,i22) %>%as_tibble() %>%filter(n_inquiry==0) %>%nrow()# 5# fix and tidy calling datacalling <-rbind(i21,i22) %>%mutate(Date=as.character(mdy(Date))) %>%mutate(group=paste(substr(Date, start=1, stop=4),Group, sep="_")) %>%# need to recreate these labels because messed up by excel in raw dataseparate(Dyad, into=c("bat_flying", "bat_roosting"), remove=F) %>%mutate(bat_flying=paste0('bat', bat_flying), bat_roosting=paste0('bat', bat_roosting)) %>%# label undirected dyadsmutate(dyad=ifelse(bat_flying<bat_roosting, paste(bat_flying, bat_roosting, sep="_"), paste(bat_roosting, bat_flying, sep="_"))) %>% dplyr::select(date= Date, group, bat_flying, bat_roosting, dyad, flight_time, n_inquiry, n_response)# combine all datad <-# start with calling data calling %>%# add SRI mutate(sri= SRI$value[match(.$dyad, SRI$dyad)]) %>%# add kinshipmutate(kinship= kinship$kinship[match(.$dyad, kinship$dyad)]) %>%# add sexes of flying and roosting batsmutate(flying.sex= sex$sex[match(.$bat_flying, sex$bat)]) %>%mutate(roosting.sex= sex$sex[match(.$bat_roosting, sex$bat)]) %>%# label sex combination for flying-->roosting dyadmutate(dyad.sex=paste(flying.sex, roosting.sex, sep="-->")) %>%# label sex combination for dyad (male, female, both)mutate(udyad.sex =case_when( flying.sex =='female'& roosting.sex =='female'~"female-female", flying.sex =='male'& roosting.sex =='male'~"male-male",TRUE~"mixed-sex")) %>%# replace zero sri (never previously seen together) with NA because association is not 0mutate(sri =if_else(sri==0, NA, sri)) %>%# if flying bat seen fewer than 4 times, then SRI is NAmutate(sri =if_else(bat_flying %in% low.n.bats, NA, sri)) %>%# if roosting bat seen fewer than 4 times, then SRI is NAmutate(sri =if_else(bat_roosting %in% low.n.bats, NA, sri))# inspect and fix cases where flying bats in more than 2 groupst <- d %>%group_by(bat_flying, date, group) %>%summarize(n=n()) %>%arrange(date) %>%group_by(bat_flying, group) %>%summarize(n=n()) %>%group_by(bat_flying) %>%summarize(n=n()) %>%filter(n>2) %>%pull(bat_flying)d %>% dplyr::select(date, group, bat_flying) %>%filter(bat_flying %in% t) %>%arrange(bat_flying) %>%as.data.frame()# some of these seem impossible and must be errors# fix impossible group labelsd$group[which(d$date=="2021-07-03"& d$bat_flying=="bat982126051278521")] <-"2021_9"d$group[which(d$date=="2021-07-03"& d$bat_flying=="bat982126058484300")] <-"2021_9"# fix cases where roosting bats in more than 2 groupst <- d %>%group_by(bat_roosting, date, group) %>%summarize(n=n()) %>%arrange(date) %>%group_by(bat_roosting, group) %>%summarize(n=n()) %>%group_by(bat_roosting) %>%summarize(n=n()) %>%filter(n>2) %>%pull(bat_roosting)d %>% dplyr::select(date, group, bat_roosting) %>%filter(bat_roosting %in% t) %>%arrange(bat_roosting) %>%as.data.frame()# fix group labelsd$group[which(d$date=="2021-07-03"& d$bat_roosting=="bat982126058484300")] <-"2021_9"# group 9 split into groups 9,1 and 9,2 during 2022# There are interesting movements between group 9 and 10 and between groups 9,1 and 9,2# For this analysis, I'm going to consider group 9,1 and 9,2 as the same groupd$group[which(d$group=="2022_9,1"| d$group=="2022_9,2")] <-"2021_9"# remove trials without inquiry or response callsd <- d %>%# 453 rowsas_tibble() %>%filter(! (is.na(n_inquiry) &is.na(n_response))) %>%# 451 rowsfilter(n_inquiry>0) %>%# 446 rowsmutate(year=substr(date, 1,4)) # add year# save data as csvwrite.csv(d, file =paste0('./processed/calling-association-kinship-data.csv'), row.names= F)```## Explore data```{r}# get dataraw <-read.csv('./data/processed/calling-association-kinship-data.csv') # look at data# head(raw)# group is the year and social group# dyad is the undirected pair (flying and roosting bat) in alphanumeric order# sri is simple ratio index of association# flying.sex is sex of flying bat# roosting.sex is sex of roosting bat# dyad.sex is the sex of the flying and roosting bat# udyad.sex is females, males, or mixed# add a few variables to the datad <- raw %>%# get log count of inquiry callsmutate(log_inquiry =log(n_inquiry)) %>%# rescale kinship so 1 unit is 0.5 mutate(kinship2 = kinship*2) %>%# if the roosting bat leaves the roost (escapes) then flight time is < 300 smutate(escape=as.numeric(flight_time<300)) %>%# convert flight time from seconds to minutes (for easier interpretation)mutate(flight_time2 = flight_time/60)###### Describe the data-------------# how many trials?# d %>% nrow()#446# # how many undirected pairs have vocal data? # d %>% # pull(dyad) %>% # unique() %>% # length#139 undirected pairs# d %>% # group_by(bat_flying, bat_roosting) %>% # summarize(n=n()) # 254 directed pairs# how many groups? # n_distinct(d$group) # 23# how many of these have association data# d %>% # group_by(dyad) %>% # summarize(sri= mean(sri)) %>% # filter(!is.na(sri))# 133 pairs have association# how many of these have kinship data# d %>% # group_by(dyad) %>% # summarize(kinship= mean(kinship)) %>% # filter(!is.na(kinship))# 82 pairs have kinship data# how many related undirected pairs# d %>% # group_by(dyad) %>% # summarize(kinship= mean(kinship)) %>% # filter(kinship>0)# 32 kin pairs# how many unrelated undirected pairs# d %>% # group_by(dyad) %>% # summarize(kinship= mean(kinship)) %>% # filter(kinship==0)# 50 nonkin pairs# what is mean kinship in group?grp.kin <- d %>%group_by(dyad, group) %>%summarize(kinship=mean(kinship, na.rm=T)) %>%group_by(group) %>%summarize(kinship=mean(kinship, na.rm=T), n=n()) %>%as.data.frame()set.seed(123)# grp.kin %>% pull(kinship) %>% boot_ci(bca=T)# 0.23, 95% CI = [0.16, 0.31]# how many groups have kin?# grp.kin %>% filter(kinship>0)# 17# how many groups have no kin?# grp.kin %>% filter(kinship==0)# 4# plot distribution of kinshipd %>%ggplot(aes(x=kinship))+facet_wrap(~udyad.sex, ncol=1, scales='free_y')+geom_histogram(fill=viridis(10, alpha =0.6)[3], color="black")+ggtitle("pairwise kinship")# # check categories# d %>%# filter(!is.na(kinship)) %>% # group_by(dyad) %>% # summarize(kinship= mean(kinship)) %>% # group_by(kinship) %>% # summarize(n=n()) ##### Effect of kinship on association-------------# plot distribution of assocd %>%ggplot(aes(x=sri))+facet_wrap(~udyad.sex, ncol=1, scale ="free_y")+geom_histogram(fill=viridis(10, alpha =0.6)[3], color="black")+xlab("association rate")+ggtitle("pairwise SRI values")# looks ok# what is mean and 95% CI of assoc among flying and roosting bats?# set.seed(123)# d %>% # group_by(dyad) %>% # summarize(sri = mean(sri)) %>% # pull(sri) %>% # boot_ci(bca=T)# 0.51, 95% CI = [0.47, 0.55]# a bit low, expected from past work is 0.76# now only nonkinset.seed(123)k1 <- d %>%filter(kinship==0) %>%group_by(dyad) %>%summarize(sri =mean(sri)) %>%pull(sri) %>%boot_ci(bca=T) %>%c(kinship=0)# k1# 0.572, 95% CI = [0.500, 0.636]# now only kin# set.seed(123)# d %>% # filter(kinship>0) %>% # group_by(dyad) %>% # summarize(sri = mean(sri)) %>% # pull(sri) %>% # boot_ci(bca=T)# 0.686, 95% CI = [0.631, 0.741]# now only close kinset.seed(123)k2 <- d %>%filter(kinship==0.5) %>%group_by(dyad) %>%summarize(sri =mean(sri)) %>%pull(sri) %>%boot_ci(bca=T)%>%c(kinship=0.5)# k2# 0.687, 95% CI = [0.628, 0.750]# now only 0.25 kinset.seed(123)k3 <- d %>%filter(kinship==0.25) %>%group_by(dyad) %>%summarize(sri =mean(sri)) %>%pull(sri) %>%boot_ci(bca=T)%>%c(kinship=0.25)# k3# 0.684, 95% CI = [0.542, 0.783]# compile meansk <-rbind(k1,k2,k3) %>%data.frame() %>%mutate(kin=as.character(kinship)) %>%mutate(kin2= kinship>0) %>%mutate(assoc= mean)# plot association by kinship(kin.dyads.plot<- d %>%filter(!is.na(kinship)) %>%group_by(dyad) %>%summarize(kinship=mean(kinship), assoc=mean(sri), udyad.sex=first(udyad.sex)) %>%mutate(kin2= kinship>0) %>%mutate(kin=as.character(kinship)) %>%ggplot(aes(x=kin, y=assoc, group= kin, color=kin2))+geom_jitter(size=2, width=0.1, height=0, aes(shape= udyad.sex))+geom_point(data= k, position =position_nudge(x =0.25), size=3)+geom_errorbar(data=k, aes(ymin=low, ymax=high, width=.1), position =position_nudge(x =0.25), size=1)+xlab("kinship")+ylab("association rate (simple ratio index)")+scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1, option ="G") +theme(legend.position="none", legend.title =element_blank()))# save as PDF# ggsave(# "kin_association.pdf",# plot = kin.dyads.plot,# scale = 1,# width = 3,# height = 5,# units = "in",# dpi = 300)##### Effect of flight time on calling ----------### How does flight time and count of inquiry calls predict count of response calls?# plot number of response calls by flight timet1 <- d %>%mutate(kinship = kinship>0.1) %>%ggplot(aes(x=flight_time, y=n_response))+geom_point(size=2, alpha=0.7, aes(color= kinship, shape= kinship))+geom_smooth(method="glm.nb", color =mako(10)[9])+scale_color_viridis_d(alpha =1, begin =0.4, end =0.9, direction =-1, option ="G") +xlab("flight time (seconds)")+ylab("response call count")+theme(legend.position="none")# plot number of inquiry calls by flight timet2 <- d %>%mutate(kinship = kinship>0.1) %>%ggplot(aes(x=flight_time, y=n_inquiry))+geom_point(size=2, alpha=0.7, aes(color= kinship, shape=kinship))+geom_smooth(method ="glm.nb", color =mako(10)[9])+scale_color_viridis_d(alpha =1, begin =0.4, end =0.9, direction =-1, option ="G") +xlab("flight time (seconds)")+ylab("inquiry call count")+theme(legend.position="none")# plot togethert <- t1+t2t# ggsave(# filename= "flight_time.pdf",# plot = t,# scale = 1,# width = 7,# height = 7,# units = c("in", "cm", "mm", "px"),# dpi = 300)```## Statistical analysis```{r}# try poisson model for effect of flight time on inquiry callingt <-glmer (n_inquiry ~ flight_time2 + (1|bat_flying) + (1|bat_roosting) + (1|group), data= d, family= poisson)check_overdispersion(t)# negative binomial model (NBM)t <-glmmTMB(n_inquiry ~ flight_time2 + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d,ziformula=~0,family=nbinom2)check_overdispersion(t)t2 <-tidy(t,conf.int=TRUE,exponentiate=T,effects="fixed", conf.method="profile")t2# for every minute of flight time, the inquiry call count increases by a factor of 1.43 [1.36, 1.49]# try poisson model for effect of flight time on response callingt <-glmer (n_response ~ flight_time2 + log_inquiry + (1|bat_flying) + (1|bat_roosting) + (1|group), data= d, family= poisson)check_overdispersion(t)# try NBMt <-glmmTMB(n_response ~ flight_time2 + log_inquiry+ (1|bat_flying) + (1|bat_roosting) + (1|group),data=d,ziformula=~1, family=nbinom2)check_overdispersion(t)check_zeroinflation(t)t2 <-tidy(t,conf.int=TRUE,exponentiate=T,effects="fixed", conf.method="profile")t2# for every minute of flight time, the response call count increases by a factor of 1.58 [1.26, 1.95]##### Effect of inquiry calling on response calling------------# plot with log countsd %>%ggplot(aes(x=log_inquiry, y=log(n_response+1)))+geom_point(size=2, color =adjustcolor("black", alpha.f =0.4))+geom_smooth(method="lm", color =mako(10)[9])+xlab("log inquiry call count")+ylab("log (x+1) response call count")# plot negative binomial curve# plot all batsd2 <- d %>%mutate(label="all pairs") %>%mutate(escape=as.logical(escape)) (t1 <-ggplot(d2, aes(x=log_inquiry, y=n_response))+facet_wrap(~label)+geom_point(size=2, alpha=0.5)+geom_smooth(method="glm.nb", color =mako(10)[9])+xlab("log of inquiry call count")+ylab("response call count"))+theme(legend.position='none')# plot with kinshipd3 <- d %>%filter(!is.na(kinship)) %>%mutate(escape=as.logical(escape)) %>%mutate(kin=if_else(kinship>0.1, "kin", "nonkin")) %>%mutate(kin=factor(kin, levels=c("nonkin", "kin"))) d2$kin <-"all pairs"shared_cols <-intersect(names(d3),names(d2))dall <-rbind(d3[, shared_cols], d2[, shared_cols])dall$kin <-factor(dall$kin, levels =c("all pairs", "nonkin", "kin")) t2 <-ggplot(dall, aes(x=log_inquiry, y=(n_response), color=kin, group=kin))+facet_wrap(~kin)+geom_point(size=2, alpha=0.5, aes(shape= escape))+geom_smooth(method="glm.nb", color =mako(10)[9])+xlab("log of inquiry call count")+ylab("response call count")+scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1, option ="G") +theme(legend.position="none")t2# save as PDF# ggsave(# "inquiry_response.pdf",# plot = inquiry.response,# width = 14,# height = 7,# units = "in",# dpi = 300)# for plots, get residuals from negative binomial mixed effect model predicting number of responses by number of inquiry callsfit <-glmmTMB(n_response ~ log_inquiry +offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d,ziformula=~1,family=nbinom2)# get difference between observed response count and predicted response count# add this response score to data for plotting# this value represents response calling more than expected from inquiry callsd$resid_response <- d$n_response -predict(fit, type="response", newdata= d, allow.new.levels = T)# plot number of inquiry calls by kinship with partnerd %>%filter(!is.na(kinship)) %>%mutate(kin=as.character(kinship)) %>%ggplot(aes(x=kin, y=n_inquiry, group= kin))+geom_violin(fill=NA, width=1)+scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1, option ="G") +geom_jitter(size=2, width=0.1, height=0, aes(color= udyad.sex))+xlab("kinship")+ylab("number of inquiry calls")+theme(legend.position="top", legend.title =element_blank())# plot number of response calls by kinship with partnerd %>%filter(!is.na(kinship)) %>%mutate(kin=as.character(kinship)) %>%ggplot(aes(x=kin, y=n_response, group= kin))+geom_violin(fill=NA, width=1)+geom_jitter(size=2, width=0.1, height=0, aes(color= udyad.sex))+scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1, option ="G") +xlab("kinship")+ylab("number of response calls")+theme(legend.position="top", legend.title =element_blank())# plot residual response call variation d %>%filter(!is.na(kinship)) %>%mutate(kin=as.character(kinship)) %>%ggplot(aes(x=kin, y=resid_response, group= kin))+geom_violin(fill=NA, width=1)+scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1, option ="G") +geom_jitter(size=2, width=0.1, height=0, aes(color= udyad.sex))+xlab("kinship")+ylab("residual response call variation")+theme(legend.position="top", legend.title =element_blank())``````{r}# plot effect of kinship on SRI(p1 <-plot_kinship_effect(d,d$sri, label='association (SRI)'))# plot effect of kinship on inquiry calling(p2 <-plot_kinship_effect(d,d$n_inquiry/(d$flight_time/60), label='inquiry calls per min of flight time'))# plot effect of kinship on response calling(p3 <-plot_kinship_effect(d,d$n_response/(d$flight_time/60), label='response calls per min of flight time'))# plot effect of kinship on response calling(p4 <-plot_kinship_effect(d,d$resid_response, label='residual response call variation'))# plot all(kinship.effects <- p1+p2+p3+p4+plot_layout(ncol =1))# save as PDF# ggsave(# "kinship_effects.pdf",# plot = kinship.effects,# scale = 1,# width = 3.5,# height = 7,# units = "in",# dpi = 300)``````{r, fig.height=10}# plot number of inquiry calls by association and kinship (by dyad type) t1 <- d %>%filter(!is.na(kinship)) %>%mutate(kin=ifelse(kinship>0.1, "kin", "nonkin")) %>%ggplot(aes(x=sri, y=log_inquiry))+facet_wrap(~dyad.sex, scales="free_y")+geom_point(aes(color=kin), size=2, alpha=0.6)+geom_smooth(method="lm", color =mako(10)[9])+xlab("association rate (simple ratio index)")+ylab("inquiry call count (log transformed)")+scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1, option ="G") +theme_bw()+theme(legend.position="none", legend.title =element_blank())# plot number of response calls by association and kinship (by dyad type)t2 <- d %>%filter(!is.na(kinship)) %>%mutate(kin=ifelse(kinship>0.1, "kin", "nonkin")) %>%ggplot(aes(x=sri, y=log(n_response+1)))+facet_wrap(~dyad.sex, scales="free_y")+geom_point(aes(color=kin), size=2, alpha=0.6)+geom_smooth(method="lm", color =mako(10)[9])+xlab("association rate (simple ratio index)")+ylab("response call count (log transformed)")+scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1, option ="G") +theme_bw()+theme(legend.position="none", legend.title =element_blank())# plot residual response calling by association (by dyad type)- controlling for inquiry callst3 <- d %>%filter(!is.na(kinship)) %>%mutate(kin=ifelse(kinship>0.1, "kin", "nonkin")) %>%ggplot(aes(x=sri, y=resid_response))+facet_wrap(~dyad.sex, scales="free_y")+geom_point(aes(color=kin), size=2, alpha=0.6)+geom_smooth(method="lm", color =mako(10)[9])+xlab("association rate (simple ratio index)")+ylab("residual response call variation")+scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1, option ="G") +theme_bw()+theme(legend.position="none", legend.title =element_blank())# combine plots(by.sex.plot <- t1+t2+t3+plot_layout(ncol =1)+plot_annotation(tag_levels ='A'))# save as PDF# ggsave(# "by_sex_plot.pdf",# plot = by.sex.plot,# scale = 1,# width = 7,# height = 14,# units = "in",# dpi = 300)```### Effect of kinship and association on response calling```{r}# poisson model is overdispersedt <-glmer(n_response ~ kinship*sri + log_inquiry +offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group), family= poisson, data = d)check_overdispersion(t)# fit negative binomial model (NBM)t <-glmmTMB(n_response ~ kinship*sri + log_inquiry +offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d,ziformula=~0,family=nbinom2)check_overdispersion(t)check_zeroinflation(t, tolerance =0.05)print("AIC")AIC(t) # 2700print("BIC")BIC(t) # 2734# fit zero-inflated NBMfit <-glmmTMB(n_response ~ kinship*sri + log_inquiry +offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d,ziformula=~1,family=nbinom2)print("AIC")AIC(fit) #2649print("BIC")BIC(fit) # 2687summary(fit)# fit zero-inflated NBM without interactionfit <-glmmTMB(n_response ~ kinship + sri + log_inquiry +offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d,ziformula=~1,family=nbinom2)print("AIC")AIC(fit) #2650print("BIC")BIC(fit) # 2684summary(fit)# get predicted valuesd$n_response_predicted <-predict(fit, type="response", newdata= d, allow.new.levels = T)# plot model performanced %>%filter(!is.na(kinship)) %>%mutate(kinship= kinship>0) %>%ggplot(aes(x=n_response, y=n_response_predicted))+geom_point(size=2, aes(color=kinship))+geom_smooth(method="lm", color =mako(10)[9])+xlab("observed count of response calls")+ylab("predicted count of response calls")+scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1, option ="G") +theme(legend.position="top")# get relationship between predicted and observedsummary(lm(n_response_predicted~n_response, data=d)) # r-squared= 0.46# get model coefficients with 95% CIscoefs <-tidy(fit,conf.int=TRUE,exponentiate=F,effects="fixed", conf.method="profile") %>%filter(component=="cond") %>% dplyr::select(-effect, -component) %>%mutate(type="full model")# coefs# save model results# write.csv(coefs, file="response_model_results.csv")# effect of SRI with kin onlyfitk <-glmmTMB(n_response ~ sri+log_inquiry +offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d %>%filter(kinship>0),ziformula=~1,family=nbinom2)coefs.k <-tidy(fitk,conf.int=TRUE,exponentiate=F,effects="fixed", conf.method="profile") %>%filter(component=="cond") %>% dplyr::select(-effect, -component) %>%mutate(type="kin only")# coefs.k# effect of SRI with nonkin onlyfitn <-glmmTMB(n_response ~ sri + log_inquiry +offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d %>%filter(kinship==0),ziformula=~1,family=nbinom2)# get 95% CIcoefs.n <-tidy(fitn,conf.int=TRUE,exponentiate=F,effects="fixed", conf.method="profile") %>%filter(component=="cond") %>% dplyr::select(-effect, -component) %>%mutate(type="nonkin only")# coefs.n# plot model resultstheme_set(theme_bw(base_size =12))plot1 <- coefs %>%filter(term !="(Intercept)") %>%mutate(term=case_when( term =='sri'~"association", term =='kinship'~"kinship", term =='log_inquiry'~"log count of inquiry calls", term =='kinship:sri'~"association x kinship interaction")) %>%mutate(term=factor(term, levels =c("log count of inquiry calls","association x kinship interaction","association", "kinship"))) %>%mutate(term =fct_rev(term)) %>%ggplot(aes(x=estimate, y=term))+geom_point(size=2)+geom_errorbarh(aes(xmin=conf.low, xmax=conf.high, height=0.2), size=1)+geom_vline(xintercept =0, linetype="dashed")+ylab("")+xlab("Effect size")+coord_cartesian(xlim=c(-1.5, 1.5))+theme(axis.text=element_text(size=12), strip.text =element_text(size=12))plot1# save as PDF# ggsave(# "response_model_results.pdf",# plot = plot1,# scale = 1,# width = 5,# height = 2,# units = "in",# dpi = 300)# check that sri and kinship do not predict response calls as single predictors # fit zero-inflated negative binomial model with only srit <-glmmTMB(n_response ~ sri + log_inquiry+offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d,ziformula=~1,family=nbinom2)summary(t)# fit zero-inflated negative binomial model with only kinshipt <-glmmTMB(n_response ~ kinship + log_inquiry+offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d,ziformula=~1,family=nbinom2)summary(t)######## Effect of kinship and association on inquiry calling---------# poisson model is overdispersedt <-glmer(n_inquiry ~ kinship*sri +log(n_response+1) +offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group), family= poisson, data = d)check_overdispersion(t)# fit negative binomial model (NBM) with log responsest <-glmmTMB(n_inquiry ~ kinship*sri +log(n_response+1) +offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d,ziformula=~0,family=nbinom2)check_overdispersion(t)AIC(t) # 3130BIC(t) # 3163# fit NBM with responsest<-glmmTMB(n_inquiry ~ kinship*sri + n_response +offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d,ziformula=~0,family=nbinom2)check_overdispersion(t)AIC(t) #3130BIC(t) # 3164summary(t)# fit NBM with responses without interactionfit2 <-glmmTMB(n_inquiry ~ kinship + sri + n_response +offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d,ziformula=~0,family=nbinom2)AIC(fit2) #3129BIC(fit2) # 3159summary(fit2)# get predicted valuesd$n_inquiry_predicted <-predict(fit2, type="response", newdata= d, allow.new.levels = T)# plot model performanced %>%filter(!is.na(kinship)) %>%mutate(kinship= kinship>0) %>%ggplot(aes(x=n_inquiry, y=n_inquiry_predicted))+geom_point(size=2, aes(color=kinship))+geom_smooth(method="lm", color =mako(10)[9])+xlab("observed count of inquiry calls")+ylab("predicted count of inquiry calls")+scale_color_viridis_d(alpha =0.8, begin =0.4, end =0.9, direction =-1, option ="G") +theme(legend.position="top")# get relationship between predicted and observedsummary(lm(n_inquiry_predicted~n_inquiry, data=d)) # r-squared= 0.71# get model coefficientssummary(fit2)coefs2 <-tidy(fit2,conf.int=TRUE,exponentiate=F,effects="fixed", conf.method="profile") %>%filter(component=="cond") %>% dplyr::select(-effect, -component) %>%mutate(type="full model")# coefs2# save model results# write.csv(coefs2, file="inquiry_model_results.csv")# plot model resultstheme_set(theme_bw(base_size =12))plot2 <- coefs2 %>%filter(term !="(Intercept)") %>%mutate(term=case_when( term =='sri'~"association", term =='kinship'~"kinship", term =='n_response'~"response calls")) %>%mutate(term=factor(term, levels =c("response calls","association", "kinship"))) %>%mutate(term =fct_rev(term)) %>%ggplot(aes(x=estimate, y=term))+geom_point(size=2)+geom_errorbarh(aes(xmin=conf.low, xmax=conf.high, height=0.2), size=1)+geom_vline(xintercept =0, linetype="dashed")+ylab("")+xlab("Effect size (log odds)")+theme(axis.text=element_text(size=12), strip.text =element_text(size=12))plot2# save as PDF# ggsave(# "inquiry_model_results.pdf",# plot = plot1,# scale = 1,# width = 7,# height = 3.5,# units = "in",# dpi = 300)# get relative amount of variance explained by random intercepttidy(fit2,conf.int=F,exponentiate=F) %>%filter(effect=="ran_pars") %>%mutate(variance= estimate^2) %>% dplyr::select(effect, group, variance) %>%group_by(effect) %>%mutate(sum=sum(variance)) %>%mutate(prop= variance/sum)# confirm no effect with single predictors# associationt <-glmmTMB(n_inquiry ~ sri + n_response +offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d,ziformula=~0,family=nbinom2)summary(t)# kinshipt <-glmmTMB(n_inquiry ~ kinship + n_response +offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),data=d,ziformula=~0,family=nbinom2)summary(t)```<!-- ' ' adds blank space --><!-- light green box --><divclass="alert alert-success"># Takeaways {.unnumbered .unlisted}- Individual ID but not kinship or association explain calling rates </div> <!-- '---' adds a gray vertical line -->--- <!-- add packages used, system details and versions --># Session information {.unnumbered .unlisted}```{r session info, echo=F}sessionInfo()```